##################################################
## Author: Stephanie L. DeMora, Jennifer L. Merolla, Brian Newman, and Elizabeth J. Zechmeister
## Project: Republican Pushback on Patriotism-Linked COVID-19 Vaccine Messages: A Note on Moral Reframing
## Purpose: Replication File for Manuscript + SI Tables & Figures; Updated probability weights
## Date: October 24th 2025
##################################################


library(haven)
library(ggplot2)
library(broom)
library(knitr)
library(kableExtra)
library(dotwhisker)
library(sjPlot)
library(sjlabelled)
library(Rmisc)
library(stringr)
library(tidyverse)
library(questionr)
library(weights)
library(jtools)
library(margins)
library(interplot)
library(janitor)
library(nnet)
library(huxtable)
library(survey)

# Data:
df <- readRDS("Replication_Data.RDS")

# Set plot theme:
steph_theme <- list(
  hrbrthemes::theme_ipsum_rc(base_size = 12,grid= "",plot_title_size = 20,axis_text_size = 18),
  ggsci::scale_color_npg(),
  ggsci::scale_fill_npg(),
  theme(legend.position = "top", plot.caption = element_text(hjust = 0, face= "italic"),
        plot.title.position = "plot",
        plot.caption.position =  "plot"))

df$treatment_2 <- relevel(df$treatment, ref = "Public Health Official")

des <- svydesign(
  ids = ~1,             # no PSUs provided
  weights = ~WEIGHT,    # general-population weight
  data = df
)

# ---------------------------------------------------------------------------------------------------------------
# ----------------------------   Manuscript Tables   ------------------------------------------------------------
# ---------------------------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------------------------
# Table 1. Vaccine and Booster Intentions
# ---------------------------------------------------------------------------------------------------------------

mod1 <- svyglm(Q4A_Likely_to_Vaccinate100 ~ treatment, design = des)
mod2 <- svyglm(Q4B_Likely_to_SecondDose100 ~ treatment, data = df, design = des)
mod3 <- svyglm(Q4C_Likely_to_Boost100 ~ treatment,  data = df, design = des)
mod4 <- svyglm(Q4D_Likely_to_Boost_OfficialRec100 ~ treatment, design = des)

export_summs(mod1, mod2, mod3, mod4,
             model.names = c("Vaccinate","Second Dose","Boost", "Boost on Recommendation"),
             stars = c(`*` = 0.1, `**` = 0.05), error_format = "(SE = {std.error}, p = {p.value})",
             note = " * p < 0.1;  ** p < 0.05",
             coefs = c('Public Health Official' = 'treatmentPublic Health Official',
                       'Democratic Voter' = 'treatmentDemocratic Voter', 
                       'Republican Voter' = 'treatmentRepublican Voter',
                       '(Intercept)' = '(Intercept)'))

# ---------------------------------------------------------------------------------------------------------------
# Table 2. Attitudes Toward COVID-19 Vaccines
# ---------------------------------------------------------------------------------------------------------------

mod5 <- svyglm(Q5A_Benefits_GreaterThan_Risks ~ treatment,  design = des)
mod6 <- svyglm(Q5B_Responsibility_to_Vaccinate ~ treatment, design = des)
mod10 <- svyglm(Q5C_Vaccines_Help_Economy ~ treatment, design = des)
mod7 <- svyglm(Q6_Family_Friend_Encourage ~ treatment, design = des)
mod8 <- svyglm(Q7_Support_Mandate_Caretakers ~ treatment, design = des)
mod9 <- svyglm(Q8_Support_Mandate_Universities ~ treatment, design = des)

export_summs(mod5, mod6, mod10, mod7, mod8, mod9,
             model.names = c("Benefits Outweigh","Responsibility","Help Economy", "Encourage", "Mandate Healthcare", "Mandate University"),
             stars = c(`*` = 0.1, `**` = 0.05), error_format = "(SE = {std.error}, p = {p.value})",
             note = " * p < 0.1;  ** p < 0.05",
             coefs = c('Public Health Official' = 'treatmentPublic Health Official',
                       'Democratic Voter' = 'treatmentDemocratic Voter', 
                       'Republican Voter' = 'treatmentRepublican Voter',
                       '(Intercept)' = '(Intercept)'))

# ---------------------------------------------------------------------------------------------------------------
# --------------------   SUPPLEMENTARY INFORMATION TABLES AND FIGURES   -----------------------------------------
# ---------------------------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------------------------
# SI Table 1.  Balance Checks
# ---------------------------------------------------------------------------------------------------------------

regression <- multinom(treatment ~ GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, data = df, na.rm = T, Hess = T, trace = F)
summary <- summary(regression)$coefficients/summary(regression)$standard.error
p <- (1 - pnorm(abs(summary), 0, 1)) * 2
p <- round(p, digits = 3)
p <- as.data.frame(p)
p$Comparison <- "Control"

#Compared to Republican Voter:
df %>%
  mutate(treatment = fct_relevel(df$treatment, "Republican Voter")) -> df_H2
regression <- multinom(treatment ~ GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, data = df_H2, na.rm = T, Hess = T, trace = F)
summary <- summary(regression)$coefficients/summary(regression)$standard.errors
p2 <- (1 - pnorm(abs(summary), 0, 1)) * 2
p2 <- round(p2, digits = 3)
p2 <- as.data.frame(p2)
p2$Comparison <- "Republican Voter"

#Compared to Democratic Voter:
df %>%
  mutate(treatment = fct_relevel(df$treatment, "Democratic Voter")) -> df_H2
regression <- multinom(treatment ~ GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, data = df_H2, na.rm = T, Hess = T, trace = F)
summary <- summary(regression)$coefficients/summary(regression)$standard.errors
p3 <- (1 - pnorm(abs(summary), 0, 1)) * 2
p3 <- round(p3, digits = 3)
p3 <- as.data.frame(p3)
p3$Comparison <- "Democratic Voter"

# Combine:
p_all <- bind_rows(p, p2, p3)
p_all %>%
  rownames_to_column("Comparison_2")

# ---------------------------------------------------------------------------------------------------------------
# SI Table 2. Descriptive Statistics by Experimental Condition
# ---------------------------------------------------------------------------------------------------------------

df %>%
  rename(Vaccinate = Q4A_Likely_to_Vaccinate100, 
         `2nd Dose` = Q4B_Likely_to_SecondDose100, 
         Boost = Q4C_Likely_to_Boost100, 
         `Boost on Recommendation` = Q4D_Likely_to_Boost_OfficialRec100,
         `Benefits Outweigh` = Q5A_Benefits_GreaterThan_Risks, 
         `Responsibility` = Q5B_Responsibility_to_Vaccinate,
         `Economy` = Q5C_Vaccines_Help_Economy,
         `Encourage` = Q6_Family_Friend_Encourage,
         `Mandate Healthcare` = Q7_Support_Mandate_Caretakers,
         `Mandate University` = Q8_Support_Mandate_Universities) %>%
  group_by(treatment) %>%
  rstatix::get_summary_stats(Vaccinate, 
                             `2nd Dose`, 
                             Boost, 
                             `Boost on Recommendation`,
                             `Benefits Outweigh`, 
                             Responsibility,
                             Economy,
                             Encourage,
                             `Mandate Healthcare`,
                             `Mandate University`,
                             type = "mean_sd") %>%
  arrange(variable)

df %>%
  tabyl(treatment)


# ---------------------------------------------------------------------------------------------------------------
# SI. Table 3. Vaccine Intentions controlling for intent to vaccinate
# ---------------------------------------------------------------------------------------------------------------

mod1 <- svyglm(Q4A_Likely_to_Vaccinate100 ~ treatment + Pre_Vax_Intentions, design = des)

huxreg('Likely Vaccinate' = mod1, 
       stars = c(`*` = 0.1, `**` = 0.05), error_format = "({std.error})", error_pos = "same",
       note = " * p < 0.1;  ** p < 0.05, standard errors in parentheses.",
       coefs = c('(Intercept)' = '(Intercept)',
                 'Democratic Voter' = 'treatmentDemocratic Voter', 
                 'Public Health Official' = 'treatmentPublic Health Official',
                 'Republican Voter' = 'treatmentRepublican Voter',
                 'Vaccination Intentions' = 'Pre_Vax_Intentions'),
       bold_signif = 0.1) 

# ---------------------------------------------------------------------------------------------------------------
# SI Table 4. Vaccine and Booster Intentions and Attitudes with Public Health Official set as baseline for comparison
# ---------------------------------------------------------------------------------------------------------------

mod1 <- svyglm(Q4A_Likely_to_Vaccinate100 ~ treatment_2, design = des)
mod2 <- svyglm(Q4B_Likely_to_SecondDose100 ~ treatment_2, design = des)
mod3 <- svyglm(Q4C_Likely_to_Boost100 ~ treatment_2, design = des)
mod4 <- svyglm(Q4D_Likely_to_Boost_OfficialRec100 ~ treatment_2, design = des)
mod5 <- svyglm(Q5A_Benefits_GreaterThan_Risks ~ treatment_2, design = des)
mod6 <- svyglm(Q5B_Responsibility_to_Vaccinate ~ treatment_2, design = des)
mod7 <- svyglm(Q6_Family_Friend_Encourage ~ treatment_2, design = des)
mod8 <- svyglm(Q7_Support_Mandate_Caretakers ~ treatment_2, design = des)
mod9 <- svyglm(Q8_Support_Mandate_Universities ~ treatment_2, design = des)
mod10 <- svyglm(Q5C_Vaccines_Help_Economy ~ treatment_2, design = des)

huxreg('Likely Vaccinate' = mod1, 
       'Likely Second Dose' = mod2, 
       'Likely Boost' = mod3, 
       'Likely Boost on Rec' = mod4,
       'Benefits Outweigh' = mod5,
       'Responsibility' = mod6,
       'Encourage' = mod7,
       'Mandate Care' = mod8,
       'Mandate Uni' = mod9,
       'Economy' = mod10,
       stars = c(`*` = 0.1, `**` = 0.05), error_format = "({std.error})", error_pos = "same",
       note = " * p < 0.1;  ** p < 0.05, standard errors in parentheses.",
       coefs = c('(Intercept)' = '(Intercept)',
                 'Democratic Voter' = 'treatment_2Democratic Voter', 
                 'Republican Voter' = 'treatment_2Republican Voter',
                 'Control' = 'treatment_2Control'),
       bold_signif = 0.1) 

# ---------------------------------------------------------------------------------------------------------------
# SI Table 5. Attitudes Toward COVID-19 Vaccines Controlling for Vaccination Status
# ---------------------------------------------------------------------------------------------------------------

mod1 <- svyglm(Q4A_Likely_to_Vaccinate100 ~ treatment + DOV_VAXSTATUS_vaxdummy, design = des)
mod2 <- svyglm(Q4B_Likely_to_SecondDose100 ~ treatment + DOV_VAXSTATUS_vaxdummy, design = des)
mod3 <- svyglm(Q4C_Likely_to_Boost100 ~ treatment + DOV_VAXSTATUS_vaxdummy, design = des)
mod4 <- svyglm(Q4D_Likely_to_Boost_OfficialRec100 ~ treatment + DOV_VAXSTATUS_vaxdummy, design = des)
mod5 <- svyglm(Q5A_Benefits_GreaterThan_Risks ~ treatment + DOV_VAXSTATUS_vaxdummy, design = des)
mod6 <- svyglm(Q5B_Responsibility_to_Vaccinate ~ treatment + DOV_VAXSTATUS_vaxdummy, design = des)
mod7 <- svyglm(Q6_Family_Friend_Encourage ~ treatment + DOV_VAXSTATUS_vaxdummy, design = des)
mod8 <- svyglm(Q7_Support_Mandate_Caretakers ~ treatment + DOV_VAXSTATUS_vaxdummy, design = des)
mod9 <- svyglm(Q8_Support_Mandate_Universities ~ treatment + DOV_VAXSTATUS_vaxdummy, design = des)
mod10 <- svyglm(Q5C_Vaccines_Help_Economy ~ treatment + DOV_VAXSTATUS_vaxdummy, design = des)


huxreg('Benefits Outweigh' = mod5,
       'Responsibility' = mod6,
       'Economy' = mod10,
       'Encourage' = mod7,
       'Mandate Care' = mod8,
       'Mandate Uni' = mod9,        
       stars = c(`+` = 0.2,`*` = 0.1, `**` = 0.05, `***` = 0.01),error_format = "({std.error})",
       coefs = c('(Intercept)' = '(Intercept)',
                 'Democratic Voter' = 'treatmentDemocratic Voter', 
                 'Public Health Official' = 'treatmentPublic Health Official',
                 'Republican Voter' = 'treatmentRepublican Voter',
                 'Vaccinated' = 'DOV_VAXSTATUS_vaxdummy'),
       bold_signif = 0.1) 


# ---------------------------------------------------------------------------------------------------------------
# SI Table 6. Vaccine and Booster Intentions and Attitudes with All Demographics as Controls
# ---------------------------------------------------------------------------------------------------------------

mod1 <- svyglm(Q4A_Likely_to_Vaccinate100 ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, design = des)
mod2 <- svyglm(Q4B_Likely_to_SecondDose100 ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, design = des)
mod3 <- svyglm(Q4C_Likely_to_Boost100 ~ treatment + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, design = des)
mod4 <- svyglm(Q4D_Likely_to_Boost_OfficialRec100 ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, design = des)
mod5 <- svyglm(Q5A_Benefits_GreaterThan_Risks ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, design = des)
mod6 <- svyglm(Q5B_Responsibility_to_Vaccinate ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, design = des)
mod7 <- svyglm(Q6_Family_Friend_Encourage ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, design = des)
mod8 <- svyglm(Q7_Support_Mandate_Caretakers ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, design = des)
mod9 <- svyglm(Q8_Support_Mandate_Universities ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, design = des)
mod10 <- svyglm(Q5C_Vaccines_Help_Economy ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, design = des)

huxreg('Likely Vaccinate' = mod1, 
       'Likely Second Dose' = mod2, 
       'Likely Boost' = mod3, 
       'Likely Boost on Rec' = mod4,
       'Benefits Outweigh' = mod5,
       'Responsibility' = mod6,
       'Encourage' = mod7,
       'Mandate Care' = mod8,
       'Mandate Uni' = mod9,
       'Economy' = mod10,
       #              stars = c(`*` = 0.1, `**` = 0.05), error_format = "({std.error})", error_pos = "same",
       #              note = " * p < 0.1;  ** p < 0.05, standard errors in parentheses.",
       coefs = c('(Intercept)' = '(Intercept)',
                 'Democratic Voter' = 'treatmentDemocratic Voter', 
                 'Public Health Official' = 'treatmentPublic Health Official',
                 'Republican Voter' = 'treatmentRepublican Voter',
                 'Gender' = 'GENDER',
                 'Age' = 'AGE4',
                 'Income' = 'INCOME4',
                 'Education' = 'EDUC5',
                 'Ideology' = 'Ideology_num'),
       bold_signif = 0.1) 


# ---------------------------------------------------------------------------------------------------------------
# SI Figure 1. Change in Vaccine Intentions and Attitudes Given Exposure to the Public Health Condition Relative to the Control Condition by Ideology (90% Confidence Intervals)
# ---------------------------------------------------------------------------------------------------------------

# Note that Interplot doesn't work with the survey package, and so we swap out
# that method for emmeans (estimated marginal means) in a more manual way below:

mod3 <- svyglm(Q4C_Likely_to_Boost100 ~ treatment * Ideology_num, design = des)
mod7 <- svyglm(Q6_Family_Friend_Encourage ~ treatment * Ideology_num, design = des)
mod8 <- svyglm(Q7_Support_Mandate_Caretakers ~ treatment * Ideology_num, design = des)
mod9 <- svyglm(Q8_Support_Mandate_Universities ~ treatment * Ideology_num, design = des)

library(emmeans)
library(dplyr)
library(purrr)
library(tidyr)
library(ggplot2)

PHO_vs_Control_curve <- function(fit, label, ideo_vals = 1:5) {
  rg <- ref_grid(fit, at = list(Ideology_num = ideo_vals))
  
  emm <- emmeans(rg, ~ treatment | Ideology_num)
  
  
  diffs <- contrast(emm, method = "trt.vs.ctrl", ref = "Control", adjust = "none")
  
  out <- summary(diffs, infer = c(TRUE, TRUE), level = 0.90) %>%
    as.data.frame() %>%
    filter(contrast == "Public Health Official - Control") %>%
    transmute(
      Ideology_num = as.numeric(Ideology_num),
      estimate,
      lower.CL,
      upper.CL,
      outcome = label
    )
  out
}

df_pho <- bind_rows(
  PHO_vs_Control_curve(mod3, "Boost"),
  PHO_vs_Control_curve(mod7, "Encourage"),
  PHO_vs_Control_curve(mod8, "Mandate Healthcare"),
  PHO_vs_Control_curve(mod9, "Mandate Universities")
)

p <- ggplot(df_pho, aes(Ideology_num, estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point() +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0) +
  facet_wrap(~ outcome, ncol = 2) +
  labs(
    x = "1 = Very liberal to 5 = Very Conservative",
    y = "Average Effect (Public Health Official – Control)",
    caption = "TESS 2022. 90% Confidence Intervals"
  ) +
  steph_theme +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 37, size = 12)
  )

p


# ---------------------------------------------------------------------------------------------------------------
# SI Fig. 2. Change in Vaccine Intentions and Attitudes Given Exposure to the Republican Voter Condition Relative to the Control Condition by Ideology (90% Confidence Intervals)
# ---------------------------------------------------------------------------------------------------------------

mod2 <- svyglm(Q4B_Likely_to_SecondDose100 ~ treatment * Ideology_num, design = des)
mod3 <- svyglm(Q4C_Likely_to_Boost100 ~ treatment * Ideology_num, design = des)
mod8 <- svyglm(Q7_Support_Mandate_Caretakers ~ treatment * Ideology_num, design = des)
mod9 <- svyglm(Q8_Support_Mandate_Universities ~ treatment * Ideology_num, design = des)

library(emmeans)
library(dplyr)
library(ggplot2)
library(purrr)
library(tidyr)

contrast_curve <- function(fit, target_label, outcome_label, ideo_vals = 1:5, conf = 0.90,
                           use_response = FALSE) {
  rg <- ref_grid(
    fit,
    at = list(Ideology_num = ideo_vals),
    type = if (use_response) "response" else "link"
  )
  
  emm <- emmeans(rg, ~ treatment | Ideology_num)
  
  diffs <- contrast(emm, method = "trt.vs.ctrl", ref = "Control", adjust = "none")
  wanted <- paste0(target_label, " - Control")
  
  summary(diffs, infer = c(TRUE, TRUE), level = conf) %>%
    as.data.frame() %>%
    filter(contrast == wanted) %>%
    transmute(
      Ideology_num = as.numeric(Ideology_num),
      estimate,
      lower.CL,
      upper.CL,
      outcome = outcome_label
    )
}

target <- "Republican Voter"

df_rv <- bind_rows(
  contrast_curve(mod2, target, "Second Dose"),
  contrast_curve(mod3, target, "Boost"),
  contrast_curve(mod8, target, "Mandate Healthcare"),
  contrast_curve(mod9, target, "Mandate Universities")
)

p_rv <- ggplot(df_rv, aes(Ideology_num, estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point() +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0) +
  facet_wrap(~ outcome, ncol = 2) +
  labs(
    x = "1 = Very liberal to 5 = Very Conservative",
    y = "Average Effect (Republican Voter – Control)",
    caption = "TESS 2022. 90% Confidence Intervals"
  ) +
  steph_theme +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 37, size = 12)
  )

p_rv


